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Abstract 

■  medial  axis  skeleton  is  a  thin  line  graph  that  preserves  the  topology  of  a  region.  The  skeleton 
often  been  cited  as  a  useful  representation  for  shape  description,  region  interpretation,  and  object 
recognition.  Unfortunately,  the  computation  of  the  skeleton  is  extremely  sensitive  to  variations  in  the 
bounding  contour.  Tiny  perturbations  in  the  contour  often  lead  to  spurious  branches  of  the  skeleton. 
In  this  paper,  we  describe  a  robust  method  for  computing  the  medial  axis  skeleton  across  a  variety  of 
scales.  The  resulting  scale-space  is  parametric  with  the  complexity  of  the  skeleton  representation.  The 
complexity  is  defined  as  the  number  of  branches  in  the  skeleton.  A  set  of  curves  is  computed  to  represent 
the  bounding  contour  across  a  variety  of  complexity  measures.  The  curves  possessing  larger  complexity 
measures  represent  greater  detail  than  urves  with  smaller  measures.  A  medial  axis  skeleton  is  computed 
directly  from  each  contour.  The  result  is  a  set  of  skeletons  that  represent  only  the  gross  structure  of  the 
region  at  coarse  scales  (low  complexity),  but  represent  more  of  the  detail  at  fine  scales  (high  complexity). 
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1  Introduction 

The  medial  axis  skeleton  is  a  thin  line  graph  that  pre¬ 
serves  the  topology  of  a  region.  The  skeleton  has  of¬ 
ten  been  cited  as  a  useful  representation  for  shape  de¬ 
scription,  region  interpretation,  and  object  recognition. 
The  skeleton  provides  a  decomposition  of  the  region  into 
salient  subparts.  It  also  provides  a  description  of  the 
connectivity  of  the  subparts. 

Unfortunately,  the  computation  of  the  skeleton  is  ex¬ 
tremely  sensitive  to  variations  in  the  bounding  contour 
of  a  region.  Tiny  perturbations  in  the  contour  often  lead 
to  spurious  branches  of  the  skeleton.  It.  is  non-trivial  to 
determine  which  of  the  branches  are  spurious  and  which 
correspond  to  significant  subregions. 

There  have  been  numerous  attempts  to  find  a  robust 
algorithm  for  computing  the  medial  axis  skeleton  (see, 
for  example,  [1],  [2],  [4],  [5],  [6],  [9],  [10],  [13]).  Most 
algorithms  use  some  deviation  of  morphological  thin¬ 
ning.  Often,  spurious  branches  are  eliminated  based 
upon  some  approximate  property  of  the  bounding  con¬ 
tour,  or  based  upon  some  property  of  the  branch  itself. 

One  common  problem  with  previous  approaches  is 
that  the  resulting  skeleton  is  inconsistent  with  the 
bounding  contour  or  the  region  from  which  it  was  com¬ 
puted.  Inconsistencies  between  the  representations  of 
the  skeleton  and  the  contour  may  lead  to  inconsistent 
inferences  in  higher  level  processes.  If  such  inconsisten¬ 
cies  could  be  eliminated,  the  performance  of  higher  level 
processes  would  be  improved. 

Another  problem  with  previous  approaches  is  that,  the 
results  are  typically  mediocre.  Most  algorithms  are  only 
capable  of  handling  simple  objects  like  pseudopods,  for 
example[9].  These  algorithms  fail  because  they  are  in¬ 
capable  of  distinguishing  between  “noise”  in  the  data 
and  subtle  features  that  may  exist  on  the  contour.  As  a 
result,  most  algorithms  tend  to  produce  skeletons  with 
spurious  branches,  or  they  tend  to  provide  skeletons  that 
are  unduly  simplified. 

Because  computation  of  the  skeleton  is  so  sensitive, 
it  is  desirable  to  represent  the  skeleton  across  a  variety 
of  scales.  The  multiple  scale  description  eliminates  the 
need  to  determine  the  “optimal”  scale  by  some  artificial 
means.  Attempts  to  find  an  optimal  scale  parameter 
in  this  and  other  contexts  typically  yield  only  marginal 
results. 

Ideally,  a  scale-space  for  the  medial  axis  skeleton 
would  provide  representations  of  the  skeleton  with  vary¬ 
ing  levels  of  detail.  At  finer  scales,  the  skeleton  would 
have  a  larger  number  of  branches;  a  greater  number  of 
features  would  be  represented.  At  more  coarse  scales,  the 
skeleton  would  have  fewer  branches;  the  skeleton  would 
represent  only  the  gross  structure  of  the  region. 

The  key  to  obtaining  a  multiple  scale  representation 
for  the  skeleton  is  to  determine  which  branches  should 
be  eliminated  as  the  algorithm  moves  from  fine  to  coarse 
scales.  Furthermore,  it  is  necessary  to  determine  the  ap¬ 
propriate  position  of  the  skeleton  branches  so  that  they 
accurately  depict  the  structure  of  the  region.  Finally,  it 
is  desirable  to  modify  the  bounding  contour  of  the  region, 
simultaneously,  so  that  each  skeleton  in  the  scale-space 
corresponds  to  a  consistent  bounding  contour. 


In  this  paper,  we  consider  a  robust  method  for  com¬ 
puting  the  medial  axis  skeleton  across  a  variety  of  scales. 

The  scale-space  is  parametric  with  the  complexity  of  the 
skeleton.  The  complexity  measure  is  defined  as  the  num¬ 
ber  of  branches  of  the  skeleton. 

The  complexity  of  the  skeleton  is  related  to  the  com¬ 
plexity  of  the  bounding  contour.  The  complexity  of  the 
contour  is  measured  by  the  number  of  extrema  of  curva¬ 
ture  contained  in  the  contour[3].  As  we  shall  see,  there  is 
a  formal  relationship  between  the  complexity  measure  of 
the  contour  and  that  of  the  skeleton.  Thus,  minimizing 
the  contour  complexity  is  tantamount  to  minimizing  the 
skeleton  complexity. 

A  set  of  curves  is  computed  to  represent  the  bound¬ 
ing  contour  across  a  variety  of  complexity  measures.  The 
curves  possessing  larger  complexity  measures  represent 
greater  detail  than  curves  with  smaller  measures.  A  me¬ 
dial  axis  skeleton  is  computed  directly  from  each  contour. 

The  result  is  a  set  of  skeletons  that  represent  only  the 
gross  structure  of  the  region  at  coarse  scales  (low  com¬ 
plexity).  but  they  represent  more  of  the  detail  at  fine 
scales  (high  complexity). 

In  Section  2,  we  discuss  the  concept  of  complexity 
in  greater  detail.  In  Section  3,  we  briefly  consider  an 
analytical  representation  for  contours;  the  contour  rep¬ 
resentation  paradigm  is  essential  for  computation  of  the 
scale-space.  In  Section  4,  we  consider  the  computation 
of  the  medial  axis  skeleton  directly  from  the  bounding 
contour.  In  Section  5,  we  define  a  scale-space  for  the 
medial  axis  skeleton  that  is  based  on  the  complexity 
measure.  In  Section  6,  we  discuss  the  benefits  of  the 
complexity  scale-space  for  the  medial  axis  skeleton  and 
compare  the  scale-space  with  the  minimum  description 
length  approach. 

2  Complexity 

The  complexity  of  an  object  may  be  viewed  as  the  num¬ 
ber  of  primitive  components  of  the  object.  Similarly,  the 
complexity  of  a  representation  of  an  object  may  be  mea¬ 
sured  by  the  number  of  subparts  contained  within  the 
representation.  In  this  section  we  seek  to  formalize  this 
notion  of  complexity  for  contours  and  the  medial  axis 
skeleton. 

Hoffman  and  Hichards[7]  have  proposed  the  use  of 
codons  to  decompose  a  contour  into  salient  parts.  They 
observe  that  minima  of  curvature  of  a  contour  serve  as 
natural  break  points  of  the  curve.  Therefore,  the  curve 
is  broken  into  sections  that  are  bounded  by  extrema  of 
curvature.  These  sections  are  called  codons.  Pairs  of 
codons  typically  correspond  to  subjective  parts  of  the 
region  bounded  by  the  contour. 

Given  this  insight,  the  number  of  codons  contained 
in  the  contour  is  a  reasonable  measure  of  the  number 
of  subjective  features  of  the  region  bounded  by  a  con-  tj 
tour.  Therefore,  the  number  of  codons  contained  in  the  Q 
contour  is  a  suitable  measure  of  the  complexity  of  the  — - 
curve.  Conveniently,  the  number  of  codons  contained  in  — - 
the  contour  is  equal  to  the  number  of  extrema  of  curva¬ 
ture  of  a  closed  contour. 

Similarly,  the  complexity  measure  of  the  medial  axis  — 
skeleton  is  the  number  of  branches  of  the  skeleton.  Each  ^3 


branch  of  the  skeleton  corresponds  to  a  subregion  of  the 
region  bounded  by  the  contour.  A  region  possessing  a 
larger  number  of  subregions  is  more  complex  than  a  re¬ 
gion  possessing  a  smaller  number  of  subregions. 

The  complexity  measure  for  contours  is  related  to  the 
complexity  measure  of  the  medial  axis  skeleton.  Each 
branch  of  the  skeleton  that  terminates  into  the  contour, 
rather  than  into  a  node  of  the  skeleton,  does  so  at  a 
positive  maximum  of  curvature.  Therefore,  an  upper 
bound  for  the  number  of  branches  in  the  skeleton  may 
be  obtained  from  the  number  of  extrema  of  curvature 
of  the  bounding  contour.  If  the  number  of  extrema  of 
a  particular  curve  is  A/,  then  there  are  at  most  M/2 
positive  maxima  of  curvature.  Therefore,  there  are  no 
more  than  M/2  branches  that  terminate  into  the  con¬ 
tour.  Each  of  these  terminal  branches  intersects  another 
branch  at  a  node  and  a  third  branch  emanates  from  the 
node.  This  branch  may  or  may  not  be  a  terminal  branch. 
In  the  worst  case,  the  number  of  non-terminal  branches 
is  M/2  —  2.  Therefore,  the  complexity  of  the  skeleton, 
B,  is  no  larger  than  M  —  2. 

A  similar  argument  may  be  made  for  a  region  with 
holes.  First,  consider  the  skeleton  associated  with  the 
bounding  contour  without  holes.  From  above,  the  skele¬ 
ton  associated  with  the  bounding  contour  has  complexity 
M  —  2.  Now  add  the  holes  one  by  one  and  consider  the 
resulting  skeleton.  Each  time  a  hole  is  added,  the  num¬ 
ber  of  branches  increases  by  no  more  than  three.  Thus, 
the  maximum  number  of  skeleton  branches  for  a  region 
with  holes  is  M  4-  3H  -2,  where  H  is  the  number  of  holes. 

More  importantly,  if  the  bounding  contour  is  de¬ 
formed  continuously  in  such  a  way  that  the  complex¬ 
ity  measure  decreases,  the  complexity  measure  of  the 
corresponding  skeleton  almost  always  decreases.  Equiv¬ 
alently,  reducing  the  number  of  extrema  of  curvature  of 
the  bounding  contour  almost  always  causes  the  number 
of  branches  of  the  skeleton  to  decrease. 

There  is  a  tradeoff  between  the  descriptive  power  of  a 
representation  and  the  complexity  of  the  representation. 
If  the  complexity  is  allowed  to  be  arbitrarily  large,  any 
set  of  data  may  be  represented.  On  the  other  hand, 
limiting  the  complexity  restricts  the  class  of  shapes  and 
objects  that  may  be  represented.  In  Section  5,  we  exploit 
this  tradeoff  to  define  a  scale-space  for  the  medial  axis 
skeleton  that  is  similar  to  a  complexity  scale-space  for 
contours.  In  the  next  two  sections, 
we  consider  an  analytical  representation  for  contours  and 
the  computation  of  the  medial  axis  skeleton  directly  from 
the  representation. 

3  Analytical  Representation  of  the 
Bounding  Contour 

The  complexity  scale-space  for  the  medial  axis  skeleton 
is  based  upon  the  complexity  scale-space  for  contours[3]. 
In  this  section,  we  briefly  consider  an  analytical  repre¬ 
sentation  for  contours  that  makes  computation  of  the 
scale-space  possible.  We  define  the  contour  represen¬ 
tation  and  consider  primitive  operations  for  deforming 
a  contour.  A  scale-space  for  contours  based  upon  the 
complexity  measure  of  the  contour  is  also  described. 


Any  reasonably  well-behaved  contour  may  be  approx¬ 
imated  by  a  list  of  pairwise  tangent  circular  arcs.  This 
representation  provides  the  ability  to  represent  the  posi¬ 
tion.  orientation,  and  curvature  of  the  contour  explicitly. 
The  representat  ion  facilitates  the  computat  ion  of  a  vari¬ 
ety  of  mathematical  properties  of  the  contour.  In  partic¬ 
ular,  the  contour  representat  ion  facilitates  t  he  computa¬ 
tion  of  an  analytical  representation  of  the  skeleton. 

We  assume  that  a  list  of  data  points  representing  the 
location  of  points  along  the  bounding  contour  is  avail¬ 
able.  The  algorithm  constructs  an  arbitrary,  initial  con¬ 
tour  that  passes  through  each  point.  The  initial  contour 
is  transformed  into  a  more  desirable  one  by  applying  a 
set  of  deformations  to  the  contour,  as  described  below. 

There  are  three  local  deformations  of  part  icular  inter¬ 
est.  The  first  operation  is  the  deformation  of  the  curva¬ 
ture  of  a  single  arc.  The  second  operation  is  the  rotat  ion 
of  two  neighboring  am.  The  third  operation  is  the  split¬ 
ting  of  a  single  arc  into  two  arcs.  These  operations  are 
illustrated  in  Figure  1. 

The  deformation  of  the  curvature  of  a  single  arc  is 
accomplished  under  the  constraint  that  its  neighboring 
arcs  remain  fixed.  This  operation  has  one  degree  of  free¬ 
dom.  As  the  radius  of  curvature  of  the  arc  changes,  the 
center  of  curvature  is  constrained  to  move  along  a  curve 
that  is  a  conic  section.  This  constraint  is  a  result  of  the 
fact  that  the  arc  of  interest  must  remain  tangent  to  its 
neighbors.  The  deformation  of  the  curvature  of  a  single 
arc  is  illustrated  in  Figure  la. 

The  rotation  of  two  neighboring  arcs  is  accomplished 
under  the  constraint  that  the  neighbors  of  the  two  rotat¬ 
ing  arcs  remain  fixed.  The  radii  of  curvature  of  the  arcs 
of  interest  also  remain  fixed;  only  the  position  of  cen¬ 
ters  of  curvature  are  modified.  Because  the  arcs  must 
remain  tangent  to  their  respective  neighbors,  the  center 
of  curvature  of  each  of  the  arcs  is  constrained  to  lie  on 
a  circle  whose  center  is  coincident  with  the  center  of  the 
respective  neighbor.  Furthermore,  the  positions  of  the 
two  arcs  of  interest  must  be  modified  in  such  a  way  that 
they  remain  tangent  to  each  other.  The  rotation  of  two 
neighboring  arcs  is  illustrated  in  Figure  lb. 

Splitting  an  arc  into  two  is  accomplished  under  the 
constraint  that  the  two  arcs  must  be  tangent  to  each 
other  and  each  of  the  arcs  is  tangent  to  one  of  the  neigh¬ 
bors  of  the  original  arc.  As  stated,  the  deformation 
has  three  degrees  of  freedom.  An  additional  constraint 
is  imposed  to  reduce  the  complexity  of  the  calculation. 
The  point  of  tangency  between  the  two  new  arcs  is  con¬ 
strained  to  lie  on  a  line  specified  by  the  algorithm.  The 
choice  of  the  constraint  line  is  dependent  on  the  context 
of  the  computation.  Splitting  an  arc  into  two  arcs  is 
illustrated  in  Figure  lc. 

The  complexity  of  a  contour,  M,  is  defined  as  the 
number  of  extrema  of  curvature  present  on  the  contour. 
Subjective  parts  of  a  region  are  delimited  by  negative 
extrema  of  curvature  on  the  bounding  contour[7].  The 
number  of  such  parts  is  limited  by  the  number  of  ex¬ 
trema  of  curvature.  Therefore,  the  number  of  extrema 
of  curvature  on  a  contour  is  related  to  the  number  of 
subjective  parts  of  the  interior  of  the  contour. 
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A  scale-space  is  constructed  based  on  the  complexity 
measure.  At  finer  scales,  the  contour  is  depicted  with  a 
higher  complexity  measure.  Thus,  more  detail  is  repre¬ 
sented.  At  more  coarse  scales,  the  contour  is  depicted 
with  a  lower  complexity  measure.  Fewer  extrema  of  cur¬ 
vature  are  present  in  the  contour.  Consequently,  less 
detail  is  present  and  only  the  gross  structure  of  the  con¬ 
tour  is  represented.  At  each  scale,  the  contour  is  chosen 
to  minimize  the  square-error  under  the  constraint  that 
it  has  the  appropriate  complexity  measure. 

The  computation  of  the  scale-space  is  performed  in 
two  stages.  In  the  first  stage,  a  curve  with  the  minimum 
complexity  measure  is  computed  under  the  constraint 
that  the  curve  passes  within  a  specified  tolerance,  £,  of 
each  data  point.  The  value,  £,  acts  as  the  scale  parame¬ 
ter.  If  <5  increases,  the  complexity  of  the  curve  decreases 
and  less  detail  is  represented.  In  the  second  stage,  the 
curve  that  is  closest  to  the  data  in  the  square-error  sense 
is  computed  under  the  constraint  that  the  curve  has  the 
complexity  measure  found  in  the  first  stage. 

In  the  first  stage,  deformations  of  the  curve  are  chosen 
to  minimize  the  difference  in  curvature  between  neigh¬ 
boring  extrema.  The  curvature  of  each  arc  associated 
with  a  maximum  is  decreased  until  it  is  not  possible  to 
do  so  without  moving  the  curve  outside  the  tolerance  of 
one  of  the  data  points.  Similarly,  the  curvature  of  each 
arc  associated  with  a  minimum  is  increased.  During  the 
course  of  this  operation,  the  number  of  extrema  is  re¬ 
duced  when  neighboring  maximum-minimum  pairs  are 
modified  such  that  their  respective  curvature  values  are 
equal. 

In  the  second  stage,  deformations  of  the  curve  are 
chosen  to  minimize  the  square-error  between  the  curve 
and  the  data  points.  The  computation  proceeds  under 
the  constraint  that  the  complexity  is  not  changed  by 
any  of  the  deformations.  During  each  iteration,  all  the 
arcs  are  modified  locally  to  reduce  the  square-error.  The 
algorithm  iterates  until  it  is  no  longer  possible  to  reduce 
the  error. 

The  result  of  the  two-stage  computation  is  the  mini¬ 
mum  complexity,  least  square-error  contour.  The  con¬ 
tour  has  the  minimum  complexity  possible  under  the 
constraint  that  the  curve  lies  within  6  of  each  data  point. 
The  curve  has  the  least  square-error  of  any  curve  that 
has  the  same  complexity  measure. 

A  multiple  scale  representation  of  the  contour  is 
achieved  by  computing  contours  with  a  variety  of  com¬ 
plexity  measures.  This  is  accomplished  by  varying  the 
scale  parameter,  6.  At  larger  values  of  6,  the  curve  has 
a  lower  complexity  measure  and  only  the  gross  structure 
of  the  contour  is  represented.  At  smaller  values  of  6,  the 
curve  has  a  higher  complexity  measure  and  more  of  the 
details  are  represented.  The  silhouette  of  an  airplane  at 
two  scales  is  depicted  in  Figure  2. 

The  analytical  representation  paradigm  provides  a  ro¬ 
bust  method  for  describing  a  contour  and  its  mathemat¬ 
ical  properties.  In  particular,  the  curvature  of  the  con¬ 
tour  is  represented  explicitly.  This  facilitates  the  compu¬ 
tation  of  a  novel  scale-space  for  contours.  Furthermore, 
the  paradigm  facilitates  the  computation  of  a  novel  scale- 
space  for  the  medial  axis  skeleton.  The  contour  represen¬ 


tation  paradigm  is  described  in  more  detail  in  an  earlier 
paper[3]. 

4  Computation  of  the  Medial  Axis 
Skeleton 

The  medial  axis  skeleton  may  be  computed  directly  from 
the  analytical  contour  representation.  In  this  section, 
we  consider  the  mechanics  of  the  computation.  First, 
we  consider  a  number  of  useful  general  properties  of  the 
skeleton  and  its  bounding  contour.  Next,  we  consider 
properties  of  the  skeleton  when  the  contour  is  made  up 
of  pairwise  tangent  circular  arcs.  Finally,  we  consider  the 
computation  of  the  skeleton  from  the  analytical  contour 
representation. 

The  medial  axis  skeleton  is  usually  defined  as  the  locus 
of  points  where  wavefronts  propagating  inward  from  the 
bounding  contour  meet  (see,  for  example,  [9] ).  The  skele¬ 
ton  points  are  locations  where  two  or  more  wavefronts 
have  propagated  the  same  distance  from  their  respective 
starting  locations.  This  definition  suggests  morpholog¬ 
ical  operators  that  approximate  the  propagation  of  the 
wavefront. 

The  medial  axis  skeleton  may  also  be  defined  by 
the  following  properties:  Each  point  on  the  skeleton  is 
equidistant  from  two  or  more  points  on  the  bounding 
contour.  There  are  no  points  on  the  boundary  closer  to 
the  skeleton  point  than  these  equidistant  points.  And, 
each  skeleton  point  lies  in  the  interior  of  the  bounding 
contour. 

This  alternate  definition  is  mathematically  equivalent 
to  the  wavefront  definition.  The  distance  from  each 
skeleton  point  to  the  closest  points  on  the  contour  is 
the  distance  traveled  by  the  associated  wavefronts.  As 
we  shall  see,  the  alternate  definition  is  constructive;  it 
leads  to  a  novel  method  of  computing  the  skeleton. 

Each  point  that  is  on  a  branch  of  a  skeleton,  but  not 
a  node,  is  equidistant  from  exactly  two  points  on  the 
contour.  Each  node  point  is  equidistant  from  three  or 
more  points  on  the  bounding  contour.  Typically,  a  node 
is  equidistant  from  exactly  three  boundary  points.  The 
case  where  the  node  is  equidistant  from  more  than  three 
points  is  a  zero  measure  condition. 

For  each  point  on  the  branch  of  the  skeleton  there  is 
a  circle  that  is  tangent  to  the  contour  in  two  places.  The 
center  of  the  circle  is  coincident  with  the  point  on  the 
branch.  The  radius  of  the  circle  is  the  distance  from  the 
center  to  the  two  nearest  points  on  the  contour.  Aside 
from  the  two  tangent  points,  the  circle  does  not  contact 
the  contour.  We  call  such  a  circle  the  interior  circle  of 
the  point  of  interest.  We  call  the  two  points  of  tangency 
between  the  interior  circle  and  the  contour  the  tangent 
points  of  the  interior  circle. 

Similarly,  for  each  node  point,  there  is  a  circle  that  is 
tangent  to  three  (or  more)  points  on  the  contour.  The 
center  of  the  circle  is  coincident  with  the  node  point  and 
the  radius  is  the  distance  from  the  node  to  the  three 
nearest  points  on  the  contour.  Aside  from  these  tangent 
points,  the  circle  does  not  contact  the  contour.  Such  a 
circle  is  called  the  interior  circle  of  the  node. 


As  a  branch  of  the  skeleton  is  traversed,  the  radius 
of  the  interior  circle  varies  continuously.  Stated  another 
way,  the  distance  from  a  skeleton  branch  to  the  contour 
varies  continuously  along  the  branch.  Furthermore,  as 
the  branch  is  traversed,  each  tangent  point  of  the  in¬ 
terior  circle  moves  continuously  along  the  contour.  In 
the  case  where  the  branch  terminates  into  the  contour, 
rather  than  into  a  node,  the  two  tangent  points  converge 
with  the  branch  at  a  point  of  maximum  curvature  on  the 
contour. 

Each  branch  of  the  skeleton  that  terminates  into  the 
contour  does  so  at  a  positive  maximum  of  curvature.  It 
is  not  the  case,  however,  that  each  positive  maximum  of 
curvature  is  associated  with  a  branch  termination.  There 
is  a  simple  test  to  determine  if  a  positive  maximum  is 
associated  with  the  termination  of  a  branch.  If  the  os¬ 
culating  circle  associated  with  the  curvature  maximum 
lies  completely  in  the  interior  of  the  contour  or  on  the 
contour,  there  is  a  branch  that  terminates  at  the  max¬ 
imum.  Otherwise,  there  is  no  terminus.  The  terminus 
test  is  illustrated  in  Figure  3. 

Now,  consider  a  contour  that  consists  of  pairwise  tan¬ 
gent  circular  arcs,  as  described  in  Section  3.  At  any  point 
on  the  skeleton,  the  tangent  points  lie  on  two  particular 
arcs  of  the  contour  representation.  Locally,  the  points  on 
the  skeleton  branch  are  equidistant  from  these  two  arcs. 
A  locus  of  points  that  is  equidistant  from  two  circles  is 
a  conic  section.  Therefore,  the  medial  axis  skeleton  con¬ 
sists  of  segments  of  curves  that  are  conic  sections.  We 
call  such  a  curve  segment  a  conic  segment  of  the  skele¬ 
ton  branch.  The  analytical  contour  representation  leads 
directly  to  an  analytic  representation  for  the  medial  axis 
skeleton. 

Because  the  conic  segments  of  the  skeleton  are  well 
characterized,  it  is  convenient  to  compute  the  skeleton  in 
a  piecewise  fashion.  The  key  to  computing  this  represen¬ 
tation  is  finding  the  end  points  of  the  branch  segments. 
At  an  end  point  of  a  branch  segment,  one  of  the  tangent 
points  of  the  interior  circle  is  guaranteed  to  be  coincident 
with  the  point  of  tangency  of  two  neighboring  contour 
arcs.  Therefore,  the  segment  end  point  must  lie  on  the 
line  determined  by  the  radius  of  the  circle  corresponding 
to  the  end  angle  of  the  arc. 

Assume  that  in  some  intermediate  stage  of  the  com¬ 
putation,  a  branch  segment  end  point  has  been  found. 
The  arcs  associated  with  the  next  branch  segment  are 
called  arcl  and  arc2,  arbitrarily.  The  line  determined 
by  the  center  of  arcl  and  the  end  point  of  arcl  is  called 
linel.  The  line  determined  by  the  center  of  arc2  and  the 
end  point  of  arc2  is  called  line2.  Assume,  without  loss  of 
generality,  that  the  branch  segment  is  hyperbolic.  There 
are  two  possibilities  for  the  location  of  the  next  segment 
end  point.  The  end  point  may  coincide  with  the  intersec¬ 
tion  of  the  hyperbola  and  linel  (candidate  pointl).  Or, 
the  end  point  may  coincide  with  the  intersection  of  the 
hyperbola  and  line2  (candidate  point2).  The  appropri¬ 
ate  choice  of  the  two  candidate  points  is  the  one  closest 
along  the  hyperbola  to  the  known  end  point.  Note  that 
the  same  reasoning  would  also  apply  to  a  branch  seg¬ 
ment  that  is  an  ellipse  or  a  parabola.  This  geometric 
situation  is  illustrated  in  Figure  4. 


The  choice  is  made  between  candidate  point  1  and  can¬ 
didate  point2  by  determining  which  point  is  closer  to 
the  previous  branch  segment  end  point  along  the  conic 
curve.  Conveniently,  there  is  a  simple  computational  test 
to  determine  the  appropriate  point.  If  candidate  pointl 
is  within  the  sector  of  arc2,  pointl  is  the  appropriate 
choice.  Similarly,  if  candidate  point2  is  within  the  sec¬ 
tor  of  arcl,  point2  is  the  appropriate  choice.  The  case 
that  both  of  these  conditions  are  true  is  zero-measure. 
Furthermore,  in  that,  case  candidate  pointl  and  candi¬ 
date  point2  are  coincident  . 

Once  the  end  points  of  the  segments  have  been  deter¬ 
mined,  it  is  possible  to  characterize  the  segment  between 
the  end  points.  The  branch  segment  is  known  to  be  a 
conic  section.  It  is  possible  to  determine  the  type  of  the 
conic  section  (hyperbola,  ellipse,  or  parabola)  by  consid¬ 
ering  the  relationship  of  the  associated  contour  arcs  and 
their  respective  curvatures.  The  foci  of  the  conic  section 
are  coincident  with  the  centers  of  the  contour  arcs.  Be¬ 
cause  the  end  points  of  the  branch  segment  lie  on  the 
conic  section,  they  provide  the  remaining  information 
necessary  to  construct  the  segment  analytically. 

Each  branch  of  the  skeleton  is  computed  in  a  piecewise 
fashion  as  described  above.  Each  segment  of  the  branch 
corresponds  to  two  arcs  on  the  curve;  the  tangent  points 
associated  with  each  point  in  the  branch  segment  lie  on 
these  two  arcs.  Two  neighboring  segments  always  share 
one  arc;  the  other  arcs  associated  with  the  two  neigh¬ 
boring  segments  are  neighbors  on  the  contour.  In  the 
example  shown  in  Figure  4,  the  segment  of  interest  is 
associated  with  arcl  and  arc2.  The  neighbor  of  this  seg¬ 
ment  is  associated  with  arcl  and  the  neighbor  of  arc2. 

In  effect,  as  the  branch  is  traversed  during  compu¬ 
tation,  the  tangent  points  on  the  contour  are  implicitly 
traversed  as  well.  The  tangent  points  associated  with 
each  point  on  a  branch  segment  are  easily  computed. 
One  of  the  tangent  points  is  simply  the  projection  of  the 
branch  point  onto  one  of  the  arcs  associated  with  the 
segment.  The  other  tangent  point  is  the  projection  of 
the  branch  point  onto  the  other  arc. 

Each  point  on  the  contour  is  associated  with  exactly 
one  point  on  the  skeleton.  It  is  not  possible  for  two 
distinct  skeleton  points  to  have  the  same  tangent  point. 
This  property  of  the  skeleton  is  useful  for  determining 
the  location  of  node  points  on  the  skeleton,  as  we  shall 
see. 

The  skeleton  computation  begins  by  determining 
starting  points  for  candidate  skeleton  branches.  Because 
it  is  known  that  branches  terminate  into  the  bound¬ 
ing  contour  at  positive  maxima  of  curvature,  these  lo¬ 
cations  are  chosen  for  the  starting  points.  As  these  can¬ 
didate  branches  are  extended,  the  locations  of  intersec¬ 
tions  of  the  branches  are  found.  At  the  intersections  of 
two  branches,  a  candidate  node  is  formed  and  an  addi¬ 
tional  candidate  branch  is  created  that  emanates  from 
the  node.  During  the  computation,  some  of  the  candi¬ 
date  branches  are  eliminated  when  it  is  determined  that 
no  branch  exists  at  its  location. 

At  each  positive  maximum  of  curvature  on  the  bound¬ 
ing  contour,  a  candidate  skeleton  branch  is  created.  By 
convention,  the  initial  branch  segment  is  the  bisecting 


radius  of  the  arc  associated  with  the  maximum  of  cur¬ 
vature.  Strictly  speaking,  such  a  segment  is  not  part 
of  the  medial  axis  skeleton  as  defined  mathematically. 
However,  these  segments  are  included  by  convent  ion  be¬ 
cause  doing  so  yields  more  intuitively  pleasing  results. 

Each  segment  is  extended  in  a  piecewise  fashion  as 
described  above.  As  the  computation  proceeds,  the  algo¬ 
rithm  must  determine  the  locations  where  the  branches 
intersect  to  form  nodes.  Each  time  a  branch  is  extended, 
the  algorithm  determines  if  the  branch  is  overextended 
relative  to  another  branch.  In  addition,  the  algorithm 
must  determine  if  the  other  branches  are  overextended 
relative  to  the  branch  of  interest.  Ultimately,  the  algo¬ 
rithm  must  determine  the  locations  of  intersections  of 
branches,  that  is,  the  nodes  of  the  skeleton.  The  follow¬ 
ing  set  of  convent  ions  achieve  these  goals. 

After  a  branch  has  been  extended  by  a  single  segment, 
the  algorithm  determines  what  interaction,  if  any,  occurs 
between  the  branch  and  the  other  branches.  Conceptu¬ 
ally,  the  algorithm  determines  if  the  t  angent  points  of  the 
branch  have  crossed  any  of  the  tangent  points  associated 
with  any  other  branch.  In  practice,  the  algorithm  con¬ 
siders  only  those  branches  that  have  arcs  of  the  contour 
in  common  with  the  branch  of  interest.  More  specifi¬ 
cally,  the  algorithm  only  considers  branches  whose  tail 
segments  (i.e.  end  segments)  have  an  arc  in  common 
with  the  tail  segment  of  the  branch  of  interest  . 

If  the  tail  segments  of  two  branches  have  a  contour 
arc  in  common,  the  systematic  application  of  a  simple 
test  determines  the  interaction  between  the  branches. 
The  test  determines  if  either  or  both  of  the  branches 
have  been  extended  beyond  the  intersection  between 
branches.  Furthermore,  these  tests  are  used  to  find  the 
node  point  which  is  located  at  the  intersection  of  two 
segments.  This  test,  described  below,  is  illustrated  in 
Figure  5. 

By  convention,  the  arc  associated  with  both  of  the 
tail  segments  is  called  the  common  arc.  We  arbitrarily 
refer  to  one  of  the  branches  as  branch  1  and  the  other  as 
branch2.  Similarly,  segment  1  and  segment‘2  are  the  cur¬ 
rent  tail  segments  of  the  respective  branches.  Arcl  and 
arc2  are  the  arcs  associated  with  segment  1  and  segment2 
that  are  not  the  common  arc. 

The  intersection  test  determines  if  an  arbitrary  point 
on  candidate  segment  1  is  beyond  the  intersection  of  seg¬ 
ment  1  and  segment2.  The  interior  circle  associated  with 
the  point  is  constructed.  By  definition  the  interior  cir¬ 
cle  is  tangent  to  arcl  and  the  common  arc;  the  center 
is  located  at  the  point  of  interest  on  segment  1.  If  the 
distance  from  the  point  to  arc2  is  greater  than  the  ra¬ 
dius  of  the  interior  circle,  the  point  of  interest  is  beyond 
the  intersection  point.  Conversely,  if  the  distance  from 
the  point  of  interest  to  arc2  is  less  than  the  radius  of  the 
interior  circle,  the  point  is  not  beyond  the  intersection  of 
segmentl  and  segments.  Of  course,  if  the  distance  from 
the  point  of  interest  to  arc2  is  equal  to  the  radius  of  the 
interior  circle,  the  point  is  the  intersection  point. 

The  intersection  test  is  illustrated  in  Figure  5.  In 
Figure  5a,  point  p?  is  beyond  the  intersection  because 
the  interior  circle  intersects  arc2.  In  Figure  5b,  point 
pi  is  not  beyond  the  intersection  because  the  interior 


circle  does  not  contact  arc2.  In  Figure  5c.  point  ;>y  is 
the  intersection  of  the  two  segments;  the  interior  circle 
is  tangent  to  all  three  contour  arcs. 

This  test  may  be  used  to  determine  if  two  tail  seg¬ 
ments  with  a  common  arc  intersect  .  The  test  is  applied 
to  both  endpoints  for  each  segment.  The  segments  inter¬ 
sect  if  and  only  if  the  intersection  test  provides  opposite 
answers  for  each  endpoint  of  both  segments.  That  is. 
one  endpoint  of  segment  1  is  beyond  the  intersection  and 
the  other  endpoint  is  not.  Similarly,  one  endpoint  of  seg- 
ment2  is  beyond  the  intersection  and  the  other  is  not. 

When  two  segments  intersect,  the  node  point  may  be 
found  using  the  intersection  test  recursively.  Initially, 
the  intersection  is  known  to  be  between  the  two  original 
endpoints  of  segmentl.  We  arbitrarily  call  these  points 
the  upper  and  lower  bound  points  of  the  node.  The  in¬ 
tersection  test  is  applied  to  a  point  midway  between  the 
bound  points.  If  the  midpoint  is  beyond  the  intersect  ion, 
the  midpoint  becomes  the  new  upper  bound.  Conversely, 
if  the  midpoint  is  not  beyond  the  intersection,  the  mid¬ 
point  becomes  the  new  lower  bound.  This  process  is 
repeated  until  bound  points  converge  to  the  node. 

It.  is  also  possible  during  the  computation  that  the 
entire  tail  segment  of  a  particular  branch  has  been  ex¬ 
tended  beyond  the  intersection  of  the  branch  ( branch  1) 
with  another  branch  (branch2).  Again,  the  intersection 
test  is  used  to  distinguish  this  situation.  The  test  is  ap¬ 
plied  to  both  endpoints  of  the  tail  segment  of  branch  1 . 
If  the  test  determines  that  both  endpoints  are  beyond 
the  intersection,  the  entire  tail  segment  is  beyond  the 
intersection.  In  that  case,  the  tail  segment  is  removed 
from  the  branch  representation. 

Note  that  extending  a  particular  branch  is  a  local  op¬ 
eration.  It  is  not  necessary  to  consider  the  entire  bound¬ 
ing  contour  to  perform  the  computation.  In  fact,  only- 
two  arcs  of  the  contour  are  required  for  each  step.  This 
suggests  that  the  branches  could  be  computed  indepen¬ 
dently,  in  parallel.  Of  course,  if  a  branch  is  extended 
such  that  its  tail  segment  has  an  arc  in  common  with 
another  tail  segment,  the  branches  must  interact  in  the 
manner  described  above.  That  is,  they  must  determine  if 
either  of  the  branches  is  overextended  and  if  the  branches 
intersect  at  a  node. 

In  our  discussion,  we  have  tacitly  assumed  that  the 
region  is  bounded  by  a  single  simply  connected  curve. 
That  is,  we  have  assumed  that  there  are  no  holes  in  the 
region.  It  is  straightforward  to  generalize  the  algorithm 
to  handle  regions  with  holes.  To  do  so  it  is  necessary 
to  find  initial  branches  such  that  each  segment  has  one 
contour  arc  on  an  interior  curve.  (An  interior  curve  is  the 
bounding  curve  of  a  hole.)  The  algorithm  extends  these 
initial  branches  to  find  nodes  similar  to  the  extension  of 
branches  described  above. 

The  initial  branches  associated  with  the  interior 
curves  are  found  in  the  following  manner:  The  point  on 
the  interior  curve  that  is  closest  to  the  bounding  curve 
is  determined.  Simultaneously,  the  point  on  the  bound¬ 
ing  contour  that  is  closest  to  the  interior  curve  is  found. 
The  midpoint  of  these  points  lies  on  the  initial  candi¬ 
date  segment  for  the  initial  branch.  An  interior  circle  is 
constructed  such  that  the  center  is  the  midpoint  and  the 


radius  is  the  distance  between  the  midpoint  and  either  of 
the  contour  points.  The  midpoint  lies  on  the  skeleton  if 
and  only  if  this  interior  circle  does  not  intersect  another 
of  the  interior  curves. 

In  the  case  that  the  interior  circle  does  intersect  an¬ 
other  interior  curve,  an  alternate  starting  point  must 
be  found.  The  alternate  starting  point  is  determined 
in  the  same  manner  described  above,  except  that  the 
closest  points  between  the  two  interior  curves  are  found 
The  midpoint  between  these  two  points  is  the  new  can¬ 
didate  skeleton  point.  The  new  candidate  point  lies  on 
the  skeleton  if  and  only  if  no  other  interior  curve  inter¬ 
sects  the  interior  circle.  In  the  event  that  the  interior 
circle  does  intersect  another  interior  curve,  the  process 
is  repeated  until  an  appropriate  skeleton  point  is  found. 
This  procedure  is  guaranteed  to  provide  a  skeleton  point 
that  is  associated  with  each  of  the  holes. 

Once  a  skeleton  point  has  been  found  for  each  hole,  a 
branch  segment  is  constructed  that,  extends  in  both  di¬ 
rections  from  each  of  the  initial  skeleton  points.  This  seg¬ 
ment  serves  as  the  starting  segment  for  an  initial  branch. 
Each  of  these  branches  is  extended  in  both  directions  to 
find  the  appropriate  nodes  with  other  branches.  Figure  6 
illustrates  the  initialization  of  a  skeleton  corresponding 
to  a  region  with  holes. 

Given  the  piecewise  conic  description  of  the  skeleton, 
it  is  possible  to  reconstruct  the  bounding  contour  ex¬ 
actly.  For  each  segment,  it  is  possible  to  reconstruct  the 
two  contour  arcs  associated  with  the  segment.  If  the 
segment  is  hyperbolic  or  elliptic,  the  centers  of  the  arcs 
are  determined  by  computing  the  locations  of  the  foci  of 
the  hyperbola  or  ellipse.  If  the  segment  is  parabolic,  the 
focus  of  the  parabola  is  the  center  point  of  one  arc;  the 
other  arc  is  a  straight  line  segment.  The  radii  of  the  two 
contour  arcs  may  be  determined  by  considering  the  radii 
of  any  two  interior  circles  along  the  conic  segment.  The 
endpoints  of  the  contour  arcs  are  determined  by  project¬ 
ing  the  endpoints  of  the  segment  onto  each  arc. 

The  medial  axis  skeleton  may  be  computed  directly 
from  the  analytical  contour  representation.  The  contour 
representation  leads  naturally  to  the  analytic  represen¬ 
tation  of  the  skeleton.  Each  branch  of  the  skeleton  is 
piecewise  elliptic,  hyperbolic,  or  parabolic.  The  bound¬ 
ing  contour  may  be  reconstructed  exactly  from  the  skele¬ 
ton. 

5  The  Medial  Axis  Skeleton 
Complexity  Scale-Space 

In  an  earlier  paper[3],  we  considered  a  complexity  scale- 
space  for  contours.  In  this  section,  we  extend  the  concept 
to  define  a  complexity  scale-space  for  the  medial  axis 
skeleton.  We  consider  the  computation  of  the  skeleton 
scale-space  from  the  contour  scale-space. 

A  scale-space  is  a  set  of  descriptions  that  differ  in 
their  level  of  detail.  At  coarse  scales  the  descriptions  are 
relatively  simple  and,  presumably,  contain  only  the  most 
important  aspects  of  the  description.  At  fine  scales,  the 
descriptions  are  relatively  complicated  and  contain  the 
details. 


A  scale-space  representation  is  desirable  because  il 
provides  alternative  descriptions  for  subsequent  process¬ 
ing.  If  a  higher  level  process  requires  accuracy  and  dense 
information,  a  fine  scale  is  appropriate.  However,  if  accu¬ 
racy  is  not  as  critical  and  sparse  information  is  sufficient, 
a  coarse  scale  is  appropriate.  In  the  latter  case,  the  com¬ 
putational  burden  is  often  significantly  reduced  because 
the  algorithm  is  required  to  process  a  smaller  quantity 
of  data. 

The  construction  of  a  scale-space  requires  a  tradeoff 
between  the  accuracy  and  the  level  of  detail  in  the  de¬ 
scription  At  fine  scales,  the  tradeoff  is  skewed  toward 
the  accuracy  of  the  desc-iption.  At  coarse  scales  the 
tradeoff  is  skewed  toward  the  simplicity  of  the  descrip¬ 
tion. 

The  measure  of  this  tradeoff  is  typically  a  smoothing 
parameter  such  as  the  spatial  width  of  a  Gaussian  filter 
applied  to  the  data  (see,  for  example,  [12]).  In  such  a 
case,  an  increased  spatial  width  of  the  filter  reduces  some 
of  the  existing  detail.  The  description  is  simplified,  but, 
the  ability  to  localize  the  remaining  components  of  the 
description  (the  accuracy)  is  reduced. 

In  this  paper,  we  propose  an  novel  method  of  quanti¬ 
fying  the  accuracy  versus  simplicity  tradeoff.  The  trade¬ 
off  yields  a  set  of  descriptions  with  varying  complexity 
measures,  as  defined  in  Section  2.  Each  description  is 
chosen  such  that  it  is  as  close  as  possible  (in  the  square- 
error  sense)  to  the  data  under  the  constraint  that  it  has 
a  particular  complexity  measure. 

In  the  case  of  a  contour,  we  assume  that  a  set  of  data 
points  along  the  contour  has  been  provided.  A  set  of 
contours  is  constructed  such  that  each  contour  has  a  dif¬ 
ferent  complexity.  That  is,  each  contour  has  a  different 
number  of  extrema  of  curvature.  Each  of  the  contours  is 
chosen  such  that  it  minimizes  the  square-error  between 
data  points  and  the  contour  under  the  constraint  that 
the  complexity  measure  is  equal  to  a  particular  value. 

In  the  case  of  the  medial  axis  skeleton,  we  also  assume 
that  a  set  of  data  points  along  the  bounding  contour  has 
been  provided.  Again,  a  set  of  skeletons  is  constructed 
such  that  the  skeletons  have  different  complexity  mea¬ 
sures.  The  bounding  contour  is  chosen  such  that  the 
square-error  between  the  data  and  the  contour  is  mini¬ 
mized  under  the  constraint  that  the  associated  skeleton 
possesses  the  appropriate  complexity  measure.  As  we 
shall  see,  the  contour  scale-space  is  very  similar  to  the 
skeleton  scale-space. 

Now,  consider  the  computation  of  the  contour  scale- 
space.  If  the  contour  is  constrained  to  pass  through  each 
data  point  exactly,  there  is  a  particular  minimum  com¬ 
plexity  measure,  Mo.  It  is  not  possible  to  construct  a 
curve  that  passes  through  every  data  point  and  has  a 
complexity  measure  smaller  than  Mo-  If  the  constraint 
is  relaxed  such  that  the  curve  must  pass  within  some  tol¬ 
erance,  b\ ,  of  each  data  poixit,  another  minimum  com¬ 
plexity  measure,  Mi,  is  obtained.  Of  course,  M\  <  Mo- 
Therefore,  as  the  tolerance,  6,  increases,  the  associated 
minimum  complexity  measure,  M ,  decreases.  Thus,  the 
tolerance,  6,  acts  as  a  scale  parameter  for  the  complexity 
scale-space. 
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For  any  tolerance,  A, ,  there  are  infinitely  many  con¬ 
tours  that  meet  the  tolerance  requirement  and  possess 
the  minimum  complexity  measure.  A/,-.  It  is  desirable 
to  choose  the  curve  that  minimizes  the  square-error  be¬ 
tween  the  data  and  the  contour  from  the  class  of  min¬ 
imum  complexity  curves.  This  suggests  a  two-stage  al¬ 
gorithm  for  determining  the  desired  minimum  complex¬ 
ity/least  square-error  curve. 

In  the  first  stage,  an  instance  of  the  minimum  com¬ 
plexity  contour  is  computed  under  the  constraint  that 
the  curve  passes  within  of  each  data  point.  The  curve 
found  in  the  first  stage  is  used  as  the  starting  point  for 
the  second  computational  stage.  The  output  of  the  first 
computational  stage  is  illustrated  in  Figure  7  with  the 
silhouette  of  an  airplane  as  test  case. 

In  the  second  comput  ational  stage,  the  minimum  com¬ 
plexity  contour  is  modified  such  that  square-error  be¬ 
tween  curve  and  the  data  points  is  minimized.  The 
modification  is  performed  under  the  constraint  that  the 
complexity  measure  does  not  change.  The  result  is  the 
minimum  complexity/least.  square-error  curve.  The  fi¬ 
nal  results  for  the  airplane  silhouette  are  depicted  in 
Figure  8.  The  computation  of  the  minimum  complex¬ 
ity/least.  square-error  curve  is  described  briefly  in  Sec¬ 
tion  3  and  more  fully  in  an  earlier  paper[3]. 

As  the  tolerance  parameter  increases  from  scale  to 
scale,  contours  with  smaller  complexity  measures  are 
found.  Some  of  the  details  in  the  contour  are  eliminated 
in  the  process;  the  more  coarse  representation  is  simpler. 
The  features  that,  are  present  in  the  representation  are 
depicted  as  accurately  as  possible.  However,  the  square- 
error  is  guaranteed  to  increase  because  the  contour  is 
unable  to  account  for  all  of  the  detail  in  the  data.  Thus, 
as  the  contour  becomes  less  complex,  the  accuracy  of  the 
representation  decreases. 

In  the  case  of  the  airplane  silhouette,  the  position  of 
the  tip  of  each  wing  is  accurately  depicted  across  the  en¬ 
tire  scale-space.  In  contrast  ,  at  coarse  scales,  the  protru¬ 
sions  on  the  back  of  the  wings  disappear  because  they 
are  smaller  than  the  scale-parameter.  The  representa¬ 
tion  is  simpler,  because  only  the  gross  structure  of  the 
wing  is  depicted.  However,  the  error  between 

the  data  and  the  contour  is  significantly  greater  in  the 
proximity  of  the  protrusions.  This  is  an  example  of  the 
tradeoff  between  simplicity  and  accuracy. 

It  would  be  reasonable  to  base  the  skeleton  scale-space 
directly  on  the  contour  scale-space.  One  option  is  to 
compute  the  contour  at  multiple  complexities.  For  each 
contour,  the  corresponding  skeleton  is  computed.  As  the 
complexity  of  the  contour  decreases,  the  complexity  of 
the  skeleton  is  guaranteed  to  decrease.  This  yields  a 
reasonable  scale-space  for  the  medial  axis  skeleton. 

However,  there  is  at  least  one  case  where  the  above 
definition  of  the  skeleton  scale-space  leads  to  a  counter¬ 
intuitive  result.  This  definition  does  not  penalize  the 
magnitude  of  the  curvature;  only  the  number  of  extrema 
of  curvature  is  considered.  Thus,  the  algorithm  prefers  a 
curve  with  relatively  large  curvature  if  such  a  choice  re¬ 
duces  the  square-error,  even  if  the  reduction  is  small.  In 
some  cases,  this  results  in  an  additional  skeleton  branch 
that  terminates  into  the  curvature  maximum.  This  phe¬ 


nomenon  is  illustrated  in  Figure  9. 

Fortunately,  a  minor  change  in  the  contour  smooth¬ 
ing  criteria  eliminates  this  effect.  The  first  stage  of  the 
contour  smoothing  algorithm  is  identical;  the  curve  is 
modified  to  minimize  the  number  of  curvature  extrema. 
In  the  second  stage,  an  additional  constraint  is  placed 
on  the  computation.  No  deformations  of  the  curve  are 
allowed  that  would  introduce  an  additional  branch  ter¬ 
minus.  The  complexity  of  the  skeleton  is  preserved  in  the 
second  stage,  as  well  as  the  complexity  of  the  contour. 

The  change  in  the  smoothing  criteria  has  only  a  small 
effect  on  the  flow  of  computation.  In  fact,  the  compu¬ 
tation  need  only  be  altered  for  arcs  corresponding  to 
positive  maxima  of  curvature.  The  computation  for  all 
other  arcs  is  identical  to  that  for  the  original  criteria. 

At  each  maximum  of  curvature  that  is  not  associated 
with  a  branch  termination,  it  is  necessary  to  ensure  that 
the  osculating  circle  extends  to  the  exterior  of  the  con¬ 
tour.  Recall  that  whenever  the  osculating  circle  does  not 
extend  to  the  exterior  of  the  contour,  a  branch  terminus 
is  guaranteed  to  occur  at  the  maximum.  Furthermore, 
whenever  the  osculating  circle  extends  to  the  exterior  of 
the  contour,  no  branch  terminus  occurs  at  the  maximum. 

The  most  obvious  method  to  avoid  the  introduction  of 
additional  skeleton  branches  is  to  disallow  deformations 
of  the  curve  that  would  create  a  spurious  branch.  How¬ 
ever,  this  method  places  an  unnecessary  burden  on  the 
computation.  It  would  be  necessary  to  test  every  can¬ 
didate  deformation  at  each  iteration  of  the  smoothing 
process  against  this  condition. 

A  more  efficient  method  is  to  compute  the  least 
square-error  curve  as  before,  then  modify  that  curve  to 
eliminate  any  spurious  branches  that  have  been  intro¬ 
duced.  This  method  requires  that  the  locations  of  the 
branch  termini  are  determined  after  the  first  stage  of  the 
smoothing  algorithm.  After  the  second  stage  of  the  con¬ 
tour  smoothing  algorithm,  each  maximum  of  curvature 
is  tested  to  determine  if  there  is  a  branch  present.  If 
a  branch  has  been  introduced  during  the  second  stage 
of  computation,  the  curve  must  be  deformed  locally  to 
eliminate  the  spurious  branch. 

To  eliminate  a  spurious  branch,  it  is  necessary  to  de¬ 
crease  the  curvature  of  the  maximum  curvature  arc.  The 
arc  is  deformed  under  the  constraint  that  the  neighbors 
remain  fixed,  as  described  in  Section  3.  The  arc  is  mod¬ 
ified  until  its  osculating  circle  contacts  another  portion 
of  the  curve;  the  branch  is  eliminated.  We  call  the  addi¬ 
tional  point  of  contact  of  the  osculating  circle  with  the 
contour  the  alternate  contact  point. 

Once  the  branch  is  eliminated,  it  is  desirable  to  mod¬ 
ify  the  curve  locally  to  reduce  the  square-error.  Such 
modifications  are  carried  out  under  the  constraint  that 
the  osculating  circle  does  not  lose  contact  with  the  con¬ 
tour  at  the  alternate  contact  point.  This  computation 
may  be  accomplished  locally  in  the  sense  that  operations 
need  information  only  about  arcs  that  are  near  the  max¬ 
imum  and  the  arc  associated  with  the  alternate  contact 
point.  No  other  arcs  of  the  contour  need  be  considered. 

The  result  of  these  computations  is  a  scale-space  rep¬ 
resentation  for  the  medial  axis  skeleton.  At  a  relatively 
coarse  scale,  a  skeleton  that  represents  only  the  gross 


structure  of  the  interior  of  the  bounding  contour  is  pro¬ 
vided.  At  a  relatively  fine  scale,  a  skeleton  that  repre¬ 
sents  more  of  the  details  is  provided.  The  scale-space  is 
illustrated  in  Figure  10  for  the  silhouette  of  an  airplane. 

In  the  case  of  the  airplane  silhouette,  the  gross  struc¬ 
ture  of  the  airplane  is  depicted  across  the  entire  scale- 
space.  At  the  coarse  scales,  the  protrusions  of  the  wings 
are  not  present  in  the  skeleton.  However,  the  gross 
structure  of  the  wings  is  accurately  represented.  At  the 
finest  scale,  the  protrusions  are  present  and  accurately 
depicted. 

Often,  the  contours  obtained  from  the  contour  scale- 
space  are  identical  with  the  bounding  contours  from  the 
skeleton  scale-space.  However,  there  are  cases  where  a 
contour  may  be  modified  such  that  an  additional  branch 
appears  (or  disappears)  in  the  skeleton  without,  changing 
the  number  of  extrema  of  curvature  in  the  curve.  In  these 
cases,  the  bounding  contour  of  the  skeleton  scale-space 
differs  slightly  from  that  of  the  contour  scale-space. 

For  each  type  of  scale-space,  a  set  of  descriptions  with 
varying  levels  of  detail  is  obtained.  The  tradeoff  between 
complexity  and  proximity  to  the  data  is  quantified  in 
terms  of  a  tolerance  about  each  data  point.  Each  rep¬ 
resentation  is  as  accurate  as  possible  in  the  square-error 
sense  under  the  constraint  of  minimum  complexity.  The 
contour  scale-space  and  the  skeleton  scale-space  are  sim¬ 
ilar;  the  difference  lies  in  a  slight  inconsistency  of  the 
complexity  measures. 

6  Discussion 

Most  “scale-space"  representations  would  be  more  accu¬ 
rately  described  by  the  term  resolution-space.  Typically, 
such  resolution-spaces  are  parameterized  by  the  spatial 
width  of  some  filter,  (usually  a  Gaussian  filter,  for  ex¬ 
ample  [12]).  Subjective  features  are  eliminated  from  the 
representation  as  the  data  is  blurred.  For  example,  at  a 
fine  scale  a  particular  feature  may  be  represented  accu¬ 
rately.  At  a  coarse  scale,  the  feature  may  be  eliminated, 
as  desired.  However,  at  an  intermediate  scale,  the  fea¬ 
ture  may  exist  with  degraded  accuracy.  There  is  no  ad¬ 
vantage  to  representing  particular  atomic  features  with 
varying  degrees  of  accuracy. 

A  true  scale-space  provides  descriptions  of  the  data 
with  varying  levels  of  detail.  At  each  scale,  all  features 
are  represented  as  accurately  as  possible.  Given  that  a 
feature  is  present  at  two  scales,  there  is  no  advantage 
in  reducing  the  accuracy  of  localization  of  the  feature 
from  one  scale  to  the  next.  It  is  more  desirable  to  retain 
the  accuracy  of  the  representation  from  one  scale  to  the 
next,  until  a  feature  is  eliminated  altogether. 

The  complexity  measure  yields  a  true  scale-space  for 
the  medial  axis  skeleton.  At  finer  scales,  greater  detail 
is  represented.  At  coarse  scales,  only  the  gross  structure 
of  the  skeleton  is  represented.  The  complexity  criterion 
explicitly  guarantees  that  the  skeleton  becomes  simpler 
at  coarse  scales. 

The  medial  axis  skeleton  complexity  scale-space  offers 
a  novel  approach  to  the  trade-off  of  accuracy  and  simplic¬ 
ity.  The  tradeoff  occurs  between  the  number  of  subjec¬ 
tive  features  and  the  square-error  between  the  data  and 
the  representation.  Whenever  a  feature  is  present  in  the 


representation,  it  is  depicted  as  accurately  as  possible. 
For  example,  the  location  of  the  tip  either  wing  of  the 
airplane  depicted  in  Figure  8  is  represented  accurately 
across  the  entire  scale-space.  Of  course,  there  is  still  a 
reduction  in  the  overall  accuracy  of  t  he  representation 
at  coarse  scales  because,  for  example,  the  protrusions  on 
the  back  side  of  each  wing  are  not  depicted  at  the  coarse 
scales. 

In  general,  it  is  more  difficult  to  obtain  a  true  scale 
space  representation  than  a  resolution-space.  It  is  non¬ 
trivial  to  obtain  a  formal  tradeoff  between  the  simplicity 
of  a  description  and  the  accuracy.  Whenever  possible, 
it  is  desirable  obtain  a  true  scale-space  representation. 
Of  course,  resolution-spaces  are  useful  to  the  extent  that 
they  approximate  the  desired  behavior  of  a  scale-space. 

Furthermore,  at  each  scale,  the  branches  that  are 
present  are  represented  as  accurately  as  possible.  That 
is,  each  branch  corresponds  to  a  portion  of  the  curve 
that  is  as  close  as  possible  to  the  data  in  the  square- 
error  sense.  The  accuracy  of  the  represent  ation  of  each 
branch  is  not  diminished  from  one  scale  to  the  next.  For 
example,  the  ability  to  localize  the  tip  of  eit  her  airplane 
wing  in  Figure  10  is  not  diminished  at  any  scale. 

Another  advantage  of  the  complexity  scale-space  for 
the  medial  axis  skeleton  is  that  the  skeleton  is  consis¬ 
tent  with  the  contour  representation.  The  mapping  be¬ 
tween  the  contour  and  the  skeleton  is  unique  and  in¬ 
vertible.  In  contrast,  most  methods  of  computing  the 
medial  axis  skeleton  yield  a  result  that  is  inconsistent 
with  the  bounding  contour;  construction  of  a  curve  from 
the  skeleton  would  lead  to  a  curve  that  differs  from  the 
original  bounding  contour. 

Consistency  among  representations  is  desirable  be¬ 
cause  higher  level  processes  may  make  inferences  based 
upon  properties  of  the  contour  or  the  skeleton.  If  the 
representations  are  consistent,  such  a  higher  level  pro¬ 
cess  is  less  likely  to  make  incompatible  inferences  about 
the  data.  Such  incompatible  inferences  would  lead  to 
degradation  of  the  overall  system  performance. 

The  complexity  scale-space  is  similar  to  the  minimum 
description  length  (MDL)  approach  at  an  intuitive  level. 
The  idea  of  providing  the  simplest  possible  representa¬ 
tion  is  present  in  both  approaches.  However,  the  for¬ 
mal  definition  of  complexity  is  different  in  the  two  ap¬ 
proaches.  As  a  result,  there  are  significant  conceptual 
differences  between  the  complexity  scale-space  and  the 
MDL  paradigm. 

The  MDL  criterion  requires  that  the  number  of  pa¬ 
rameters  employed  by  a  model  to  account  for  the  data 
should  be  minimized[llj.  More  precisely,  the  number  of 
bits  required  to  encode  the  data  is  minimized.  Thus, 
an  important  component  of  the  MDL  approach  is  the 
choice  of  an  efficient  model  for  the  representation.  For¬ 
mally,  the  MDL  approach  requires  that  the  represen¬ 
tation  correspond  to  an  optimal  code  in  the  information 
theoretic  sense.  A  theory  for  choosing  the  representation 
from  a  priori  probability  distributions  is  well  known  (see 
Leclerc[8],  for  example). 

Unfortunately,  the  determination  of  such  an  optimal 
code  for  a  general  application  is  difficult  in  practice. 
Typically,  minimum  description  length  approaches  as- 


sume  a  particular  form  of  the  representation.  The  num¬ 
ber  of  parameters  employed  by  the  representation  is  min¬ 
imized.  Interestingly,  the  MDL  theory  provides  a  simple 
objective  measure  of  the  performance  of  a  particular  rep¬ 
resentation.  The  performance  measure  is,  of  coarse,  the 
average  length  of  the  code. 

For  example,  Leclerc[8]  applies  the  MDL  technique  to 
the  image  partitioning  problem.  Each  subregion  requires 
a  set  of  parameters  to  distinguish  the  boundary  and  to 
specify  the  behavior  of  image  brightness  within  the  sub- 
region.  The  additional  parameters  required  for  each  sub- 
region  causes  the  algorithm  to  favor  a  small  number  of 
partitions.  However,  when  there  are  relatively  few  par¬ 
titions,  it  is  more  difficult  to  account  for  the  variations 
from  the  model  within  each  subregions;  a  larger  number 
of  bits  per  subregion  is  required.  Conversely,  this  effect, 
favors  a  larger  number  of  partitions.  There  is  a  partic¬ 
ular  choice  of  the  partitioning  that  optimizes  these  two 
opposing  effects. 

In  contrast,  the  complexity  scale-space  seeks  to  mini¬ 
mize  the  number  of  subjective  features  in  the  description. 
The  optimization  is  indifferent  to  the  number  parame¬ 
ters  used  by  the  representation.  In  the  case  of  a  con¬ 
tour,  the  contour  may  be  described  by  arbitrarily  many 
circular  arcs;  the  complexity  depends  only  on  the  num¬ 
ber  of  extrema  of  curvature.  Similarly,  each  branch  of 
a  medial  axis  skeleton  may  be  represented  by  arbitrar¬ 
ily  many  branch  segments;  the  complexity  is  specified 
by  the  number  of  branches.  An  MDL  approach  would 
essentially  require  minimizing  the  number  of  arcs  in  a 
contour  or  the  number  of  segments  in  a  skeleton. 

Furthermore,  the  complexity  scale-space  does  not  at¬ 
tempt  to  provide  an  efficient  code  for  the  representation. 
Rather,  the  scale-space  seeks  to  make  the  most  relevant 
information  explicit  for  higher  level  processes.  In  fact, 
there  is  a  loss  of  information  at  the  coarse  scales.  We 
have  argued,  above,  that  this  loss  of  information  is  de¬ 
sirable  for  computational  purposes. 

Another  important  distinction  is  that  MDL  ap¬ 
proaches  typically  do  not  provide  scale-space  representa¬ 
tions.  MDL  approaches  provide  a  single  representation 
for  each  data  set.  In  effect,  an  optimal  scale  is  chosen 
implicitly  by  the  MDL  criterion.  The  scale  chosen  by 
the  MDL  criterion  is  the  one  requires  the  fewest  number 
of  bits  to  represent  the  data. 

The  paradigm  described  in  this  paper  leads  to  true 
scale-spaces  for  contours  and  the  medial  axis  skeleton. 
At  first  glance,  the  complexity  scale-spaces  are  very  sim¬ 
ilar  to  minimum  description  length  approaches;  the  no¬ 
tion  of  simplicity  of  the  representation  is  exploited  in 
each  case.  However,  there  are  significant  philosophical 
and  practical  differences  between  the  complexity  scale- 
spaces  and  MDL  approaches. 

7  Conclusion 

We  have  considered  a  novel  approach  to  computing  the 
medial  axis  skeleton  across  a  variety  of  scales.  Complex¬ 
ity,  as  defined  by  the  number  of  branches  of  the  skeleton, 
is  a  natural  measure  of  the  level  of  detail  of  the  skeleton 
description.  The  multiple  complexity  paradigm  leads  to 
a  true  scale-space;  the  ability  to  localize  a  feature  is  not 


diminished  from  one  scale  to  another  unless  t  he  feature  is 
completely  eliminated.  The  representation  for  the  skele¬ 
ton  is  consistent  with  the  corresponding  representation 
of  the  bounding  contour;  the  mapping  between  the  skele¬ 
ton  and  the  bounding  contour  is  unique  and  invertible. 

References 

[1]  D.  H.  Ballard  and  C.  M.  Brown.  Computer  Vision 
Prentice- Hall,  Englewood  Cliffs.  NJ.  1982. 

[2]  H.  Blum.  A  transformation  for  extracting  new  de¬ 
scriptions  of  shape.  In  Symposium  on  Modi  Is  for  the 
Perception  of  Speech  and  Visual  Form.  MIT  Press, 
Cambridge,  MA,  1964. 

[3]  R.  D.  Chaney.  Analytical  representation  of  con¬ 
tours.  AI  Memo  1392,  MIT  Artificial  Intelligence 
Laboratory.  October  1992. 

[4]  R.  T.  Chin,  H.  K.  Wan,  D.  L.  Stover,  and  R.  D.  lver- 
sion.  A  one-pass  thinning  algorithm  and  its  parallel 
implementation.  Computer  Vision.  Graphics,  and 
Image  Processing.  40(  1  );30— 40,  1987. 

[5]  A.  R.  Dill,  M.  D.  Levine,  and  P.  B.  Nobel.  Multiple 
resolution  skeletons.  Computer  Vision.  Graphics, 
and  Image  Processing ,  40(  1  ):30— 40,  1987. 

[6]  S.  B.  Ho  and  C.  R.  Dyer.  Comparison  of  thinning 
algorithms  on  a  parallel  processor.  IEEE  Transac¬ 
tions  on  Pattern  Analysis  and  Machine  Intelligence, 
8(4):512-520,  1986. 

[7]  D.  D.  Hoffman  and  W.  A.  Richards.  Parts  of  recog¬ 
nition.  Cognition,  pages  65-96,  1984.  Look  this  up. 

[8]  Y.  G.  Leclerc.  Constructing  simple  stable  descrip¬ 
tions  for  image  partitioning.  International  Journal 
of  Computer  Vision,  3(  1):73— 102,  May  1989. 

[9]  F.  Leymarie  and  M.  D.  Levine.  Simulating  the 
grassfire  transform  using  an  active  contour  model. 
IEEE  Transactions  on  Pattern  Analysis  and  Ma¬ 
chine  Intelligence,  14(  1  ):56— 75,  Jan.  1992. 

[10]  M.  P.  Martinez-Perez,  J.  Jimenez,  and  J.  Navalon. 
A  thining  algorithm  based  on  contours.  Computer 
Vision,  Graphics,  and  Image  Processing,  39(2):186- 
201,  August  1987. 

[11]  J.  Rissan.  Minimum-description-length  principle. 
In  Encyclopedia  of  Statistical  Sciences,  volume  5, 
pages  523-527.  Wiley,  New  York,  1987. 

[12]  A.  P.  Witkin.  Scale  space  filtering.  In  Proc.  7th 
International  Conference  on  Artificial  Intelligence, 
pages  1019-1021,  Karlsruhe,  1983. 

[13]  Y.  Xia.  Skeletonization  via  the  realization  of  the  fire 
front's  propagation  and  extinction  in  binary  shapes. 
IEEE  Transactions  on  Pattern  Analysts  and  Ma¬ 
chine  Intelligence,  11(10);  1076-1086,  October  1989. 


9 


Figure  1:  Local  deformations  of  an  analytic  contour.  In  (a),  the  curvature  of  a  single  arc  is  modified  under  the 
constraint  that  it  remains  tangent  to  its  neighbors.  In  (b),  two  neighboring  arcs  are  rotated  under  the  constraint 
that  they  remain  tangent  to  their  neighbors  and  to  each  other.  In  (c),  an  arc  is  split  into  two  arcs.  The  two  new 
arcs  are  tangent  to  each  other  at  a  point  along  a  specified  constraint  line. 


Figure  2:  Contours  at  two  different  scales.  In  (a)  a  fine  scale  representation  of  the  silhouette  of  an  airplane  is 
depicted,  in  (b)  a  coarse  scale  representation.  The  small  squares  near  each  contour  represent  the  original  data  points. 
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(a)  (b) 


Figure  3:  Test  for  the  existence  of  branch  terminating  into  a  maximum  of  curvature.  Each  illustration  depicts  a 
portion  of  a  contour  and  its  medial  axis  skeleton;  the  osculating  circle  for  the  maximum  of  curvature  at  the  top  of  the 
curve  is  also  shown.  In  (a),  no  branch  terminates  into  maximum  of  curvature  because  the  osculating  circle  extends 
to  the  exterior  of  the  curve.  In  (b),  a  branch  does  terminate  into  the  maximum  of  curvature  because  the  osculating 
circle  remains  in  the  interior  of  the  contour. 


Figure  4:  Computation  of  a  branch  segment.  The  branch  segment  is  the  hyperbolic  curve  connecting  point  pO  to 
point  p2.  The  hyperbola  is  defined  by  the  property  that  each  point  is  equidistant  from  arcl  and  arc2  (two  arcs  in  the 
contour  representation).  The  candidate  end  point,  pi,  is  the  intersection  of  the  end  radius  of  arcl  and  the  hyperbola. 
Similarly,  the  candidate  end  point,  p2,  is  the  intersection  of  the  end  radius  of  arcl  and  the  hyperbola.  The  point, 
p2,  is  chosen  because  it  is  within  the  sector  of  arcl  -  the  point,  pi,  is  not  in  the  sector  of  arc2.  The  next  segment 
that  would  be  computed  extends  beyond  point,  p2,  and  is  determined  by  arcl  and  the  neighbor  of  arc2. 
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Figure  5:  Test  for  determining  the  location  of  node  points.  In  each  figure,  three  segments  that  intersect  at  a  node 
are  depicted.  The  contour  arcs  associated  with  these  segments  are  also  depicted.  Segmentl  is  equidistant  from  the 
common  arc  and  arcl;  segment2  is  equidistant  from  th  common  arc  and  arc2;  and  segment3  is  equidistant  from 
arcl  and  arc2.  Conceptually,  segmentl  and  segment2  intersect  to  form  a  node;  segment3  emanates  from  this  node. 
In  (a)  point  p2  is  beyond  the  intersection  of  segmentl  and  segment2:  the  interior  circle  is  tangent  to  the  common 
arc  and  the  interior  circle  intersects  arc2.  In  (b),  point  p\  on  segmentl:  the  interior  circle  is  tangent  to  the  common 
arc  and  arcl,  and  the  interior  circle  does  not  intersect  arc2.  In  (c),  point  pz  is  the  intersection  between  segmentl 
and  segment2:  the  interior  circle  is  tangent  to  all  three  contour  arcs. 
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Figure  6:  Skeleton  for  a  region  with  holes.  The  computation  of  a  skeleton  for  a  region  with  two  holes  is  shown.  In  (a), 
the  initial  skeleton  points  associated  with  the  interior  curves  are  depicted  along  with  their  respective  interior  circles. 
In  (b),  the  initial  segments  are  extended  from  the  initial  points.  These  segments  serve  to  initialize  the  skeleton 
branches  for  the  computation.  In  (c),  the  skeleton  associated  with  the  region  is  shown. 
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Figure  9:  Sensitivity  of  the  skeleton  at  maxima  of  curvature.  Two  curves  are  shown  with  their  respective  medial 
axis  skeletons.  The  two  curves  are  identical  except  near  the  maximum  of  curvature  on  the  top  of  each  curve.  The 
complexity  measures  of  the  curves  are  equal.  The  top  curve  has  an  additional  skeleton  branch  because  the  magnitude 
of  curvature  at  the  maximum  is  much  greater  than  that  of  the  bottom  curve.  The  complexity  measure  of  the  skeleton 
is  very  sensitive  to  subtle  changes  in  the  bounding  contour,  particularly  near  positive  maxima  of  curvature. 
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Figure  10:  The  complexity  scale-space  for  the  medial  axis  skeleton.  The  medial  axis  skeleton  and  the  corresponding 
bounding  contour  are  shown  at  four  scales.  The  scale  parameter,  6,  is  doubled  between  each  scale.  Note  that  the 
number  of  branches,  as  well  as  the  number  of  features  on  the  contour,  decreases  as  the  scale  parameter  increases. 
The  features  that  are  present  in  each  representation  are  depicted  as  accurately  as  possible  in  the  square-error  sense. 
The  small  squares  near  each  contour  represent  the  initial  data  points  of  the  contour. 
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